home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / program / nrpas13.zip / SVDCMP.DEM < prev    next >
Text File  |  1991-04-29  |  2KB  |  95 lines
  1. PROGRAM d2r9 (input,output,dfile);
  2. (* driver for routine SVDCMP *)
  3. LABEL 10,99;
  4. CONST
  5.    np = 20;
  6.    mp = 20;
  7. TYPE
  8.    glnparray = ARRAY [1..np] OF real;
  9.    glmpbynp = ARRAY [1..mp,1..np] OF real;
  10.    glnpbynp = ARRAY [1..np,1..np] OF real;
  11. VAR
  12.    j,k,l,m,n : integer;
  13.    a,u : glmpbynp;
  14.    v : glnpbynp;
  15.    w : glnparray;
  16.    dfile : text;
  17.  
  18. (*$I MODFILE.PAS *)
  19. (*$I SVDCMP.PAS *)
  20.  
  21. BEGIN
  22. (* read input matrices *)
  23.    glopen(dfile,'matrx3.dat');
  24. 10:   readln(dfile);
  25.    readln(dfile);
  26.    readln(dfile,m,n);
  27.    readln(dfile);
  28. (* copy original matrix into u *)
  29.    FOR k := 1 to m DO BEGIN
  30.       FOR l := 1 to n DO BEGIN
  31.          read(dfile,a[k,l]);
  32.          u[k,l] := a[k,l]
  33.       END;
  34.       readln(dfile)
  35.    END;
  36.    IF (n > m) THEN BEGIN
  37.       FOR k := m+1 to n DO BEGIN
  38.          FOR l := 1 to n DO BEGIN
  39.             a[k,l] := 0.0;
  40.             u[k,l] := 0.0
  41.          END
  42.       END;
  43.       m := n
  44.    END;
  45. (* perform decomposition *)
  46.    svdcmp(u,m,n,np,np,w,v);
  47. (* write results *)
  48.    writeln ('Decomposition matrices:');
  49.    writeln ('Matrix u');
  50.    FOR k := 1 to m DO BEGIN
  51.       FOR l := 1 to n DO BEGIN
  52.          write (u[k,l]:12:6);
  53.       END;
  54.       writeln
  55.    END;
  56.    writeln ('Diagonal of matrix w');
  57.    FOR k := 1 to n DO BEGIN
  58.       write(w[k]:12:6)
  59.    END;
  60.    writeln;
  61.    writeln ('Matrix v-transpose');
  62.    FOR k := 1 to n DO BEGIN
  63.       FOR l := 1 to n DO BEGIN
  64.          write(v[l,k]:12:6)
  65.       END;
  66.       writeln
  67.    END;
  68.    writeln;
  69.    writeln ('Check product against original matrix:');
  70.    writeln ('Original matrix:');
  71.    FOR k := 1 to m DO BEGIN
  72.       FOR l := 1 to n DO BEGIN
  73.          write (a[k,l]:12:6)
  74.       END;
  75.       writeln
  76.    END;
  77.    writeln ('Product u*w*(v-transpose):');
  78.    FOR k := 1 to m DO BEGIN
  79.       FOR l := 1 to n DO BEGIN
  80.          a[k,l] := 0.0;
  81.          FOR j := 1 to n DO BEGIN
  82.             a[k,l] := a[k,l]+u[k,j]*w[j]*v[l,j]
  83.          END
  84.       END;
  85.       FOR l := 1 to n-1 DO write(a[k,l]:12:6);
  86.       writeln(a[k,n]:12:6);
  87.    END;
  88.    writeln ('***********************************');
  89.    IF eof(dfile) THEN GOTO 99;
  90.    writeln ('press RETURN for next problem');
  91.    readln;
  92.    GOTO 10;
  93. 99:   close(dfile)
  94. END.
  95.